† Corresponding author. E-mail:
Project supported by the National Natural Science Foundation of China (Grant Nos. 11574310, 11674345, and 21733010) and Beijing National Laboratory for Molecular Sciences, China (Grant No. BNLMS201835).
It is very important to determine the phase transition temperature, such as the water/ice coexistence temperature in various water models, via molecular simulations. We show that a single individual direct simulation is sufficient to get the temperature with high accuracy and small computational cost based on the generalized canonical ensemble (GCE). Lennard–Jones fluids, the atomic water models, such as TIP4P/2005, TIP4P/ICE, and the mW water models are applied to illustrate the method. We start from the coexistent system of the two phases with a plane interface, then equilibrate the system under the GCE, which can stabilize the coexistence of the phases, to directly derive the phase transition temperature without sensitive dependence on the applied parameters of the GCE and the size of the simulation systems. The obtained result is in excellent agreement with that in literatures. These features make the GCE approach in determining the phase transition temperature of systems be robust, easy to use, and particularly good at working on computationally expensive systems.
The first-order phase transitions, such as the water/ice transition, are very popular phenomena in nature. Molecular modeling and simulations offer microscopic insights to both the thermodynamics and kinetics of phase transitions.[1–11] The phase coexistence temperature is one of the center values to verify models, and efficiently calculating the temperature via simulations to compare with experiments is key in the molecular modeling and simulations of phase transitions.
Estimating the phase transition temperature via simulations has been a long term issue. Take the TIP4P water model for example, it was proposed by Jorgensen in 1983[12] and well performs in revealing the density of liquid water, the density near the critical point, and the enthalpy of vaporization.[13] A rough estimate of the melting point of ice Ih was first given by Kroes[14] as 230 K < Tc < 250 K for the water model by monitoring the translational mobility of atoms on the ice/water interfaces. Later in 2000, with free energy (FE) calculation, Gao indicated Tc = 238 ±7 K.[6] While another study proposed by Vega in 2005 implied a similar result, Tc = 232 ±5 K.[15] Although the FE approach determined phase transition temperature is theoretically flawless, refers to the state at which two phases have Gibbs free energy identical, tiny uncertainty in free energy estimation can lead to a totally different result. Gao described this as “1 percent in relative error leads to a shift of melting point more than 10 K”, the FE approach seems to reach its limit.
Direct coexistence (DC) simulation of phase transition is also a generally accepted method to determine the phase transition temperature. This method is particularly popular in extracting the gas–liquid coexisting line on the phase diagram,[16–19] and also can be applied in liquid/solid systems.[5,20,21] Those works of determining the phase transition temperature with DC approach show great agreement with the FE approach, and sometimes better precision if the data sample is sufficient. For example, with the water model TIP4P-2005 (by reparameterising from the TIP4P),[22] Conde calculated the melting point of the ice Ih by both FE and DC approaches, that is, 252±5 K and 249±3 K, respectively.[21] In 2017, by performing simulations covering a range of temperature space, repeating the simulation with different seeds for 5 times, using the potential energy evolution in tens of nanoseconds as the phase transition indicating order parameter, Conde developed the DC approach with exhaustive study, obtaining Tc = 249.5 ± 0.1 with the best precision to date.[5] However, the plenty of simulations and expensive CPU cost have limited the application of this approach.
The Hamiltonian Gibbs–Duhem integration offers another option. A number of studies have demonstrated that once the melting point of a model is known, the melting point of a different model can be easily estimated with the Gibbs–Duhem integration.[10,15] However, this reference model dependent estimate is not universally applicable, one can neither extract the melting point of ice Ic from an ice Ih reference, nor a silicon phase transition reference.
In this work, we aim to provide a simple and universal method with high precision and efficiency for determining the phase transition temperature. To achieve this, the generalized canonical ensemble (GCE)[23–28] has been implemented successfully in Lennard–Jones (LJ) particles and ice/water systems with the mW, TIP4P-2005, TIP4P-ICE water models. Unlike in the normal canonical ensemble simulations where phase-coexistence states are unstable thus can not be focused on, the GCE allows sufficiently sampling in the phase-coexistence regions, making a single simulation of DC approach be enough to calculate the phase transition temperature with considerable precision.
This work is organized as follows. Section
The GCE method has been explicitly demonstrated in previous studies,[23–25] here we present a brief overview of this method. In GCE, the energy-dependent thermostat temperature is applied to couple with the system, in physics, it is equivalent to apply a finite-size thermostat so that the temperature of the thermostat varies as the energy of the system (and the energy of the thermostat, since the total energy of the system and thermostat is constant). In GCE, we choose the simplest form of the energy-dependent (reverse) temperature
![](cpb_29_8_080505/cpb_29_8_080505_ieqn1.gif)
The GCE (or gNPT) is equivalent to applying an effective potential energy, thus its implementation is simple and direct by rescaling the physical force on each atom. It is also possible to reset the temperature of the thermostat at each MD simulation step under the original potential energy surface, from the viewpoint of applying a finite-size thermostat. The two implementations are thought to be equivalent to each other for achieving the equilibrium properties.[25,26]
The GCE aims to use a large positive α to visit sufficiently the phase-coexistent conformational regions, while the canonical ensemble (α = 0) can not focus in the regions due to the thermodynamical instability of the coexistent states. We define the microcanonical entropy S(E) = ln ∫ drδ(E – U(r)), and will show that S″(E) = ∂2 S/∂ E2 > 0 is the condition of the thermodynamical instability. The probability distribution of GCE in the energy space is
If neglecting the possible small asymmetry of the GCE probability distribution around Em, we approximately have
At the phase transition temperature (or named as the phase coexistence temperature) Tc, the two pure phases have the same Helmholtz (or Gibbs) free energy as a function of temperature, As(Tc) = Al(Tc). Here we use the Landau free energy F(E; T) = E – T S(E), and
![](cpb_29_8_080505/cpb_29_8_080505_ieqn8.gif)
![](cpb_29_8_080505/cpb_29_8_080505_ieqn9.gif)
![](cpb_29_8_080505/cpb_29_8_080505_ieqn10.gif)
Usually, ones can run simulations in canonical ensembles to estimate the temperature Tm at which the system often melts from solid to liquid phase, and the temperature Tf at which the system easily freezes from liquid to solid. The coexistence temperature Tc is between the two temperatures, and is possibly estimated as
In GCE, the whole S′(E) curve can be constructed, thus we can integrate to get S(E) and to estimate Tc from Eq. (
We assess the method in determining the coexistence temperature of the Lennard–Jones fcc crystal, and that of ice Ih in the all-atom TIP4P-2005 and TIP4P-ICE water models, as well as the coarse-grained mW water.[7]
The simulation system is built with 32480 or 4584 water molecules for mW (approximately 400 Å × 56 Å × 43 Å) and all-atom water models (approximately 200 Å × 26 Å × 29 Å), respectively. The ice phase is generated from 20 ns equilibrated bulk ice Ih at 260 K, with x–y plane set as ice Ih basal plane, leaving the secondary prismatic plane (
All MD simulations are carried out in the GCE embedded isothermal-isobaric (NPT) ensembles (gNPT), employing LAMMPS molecular dynamics code with our added GCE module.[23] A Nose–Hoover barostat[29] with 1 ps relaxation time is used to keep the pressure at 1 atm. For all-atom water models, a time step of 2 fs is used in all simulations. The pair interaction between the oxygen atoms is truncated smoothly at 13.0 Å. Nontruncated electrostatic interactions are treated by the particle–particle particle mesh solver (pppm) with a real space cutoff of 13.0 Å and precision tolerance of 10−5. The simulation parameters in the mW water model are set normally, similar to that in reference.[7] In the LJ systems, a cutoff of 2.5σ is used to truncate the interaction.
We first illustrate the gNPT simulations of the coexistent solid/liquid LJ system to determine the phase transition temperature and its dependence on the parameters of gNPT.
The study starts with a set of gNPT simulations of LJ particles with varied parameter H0/N to visit the states with different regions of enthalpy, corresponding to the different temperatures in pure phases and the coexistent phases with different fraction of two phases. The obtained mean enthalpy of each particle in every gNPT simulation, H/N, is from –7 to –5 (unit is the ε of the LJ potential). Here N is the number of particles. The result can be easily identified as a pure phase or a phase-coexistent state, as shown in Fig.
The coexistent temperature Tc is approximately extracted as the temperature in the middle of the coexistent region, or as the statistical average of the corresponding temperature from the platform region (–6.4 ⩽ H/N ⩽ –5.8).
We also carry out the gNPT simulations with water models. Water/ice coexisting systems are first simulated at 250 K for 100 ps to relax the interfaces and to estimate the parameter of gNPT, H0/N. For the mW model, the mean enthalpy in the coexistent region is around –10.7 (the unit is always kcal/mol in this work for water). Then we run gNPT simulations each for 20 ns with the mean enthalpy ranged from –9.9 to –11.5 are shown in Fig.
It is interesting to see the temperature evolution in the gNPT simulations, as shown in Fig.
![]() | Fig. 4. The temperatures obtained in the gNPT simulations. Cyan region refers to 274.6 ± 1 K, the melting point of the mW water model proposed in previous study.[7] Tc = 275.1 ± 0.3 K is the average of six simulations located in the phase-coexistence region. |
It has to be point out that the melting temperatures we extracted from the gNPT ensembles at different assigned enthalpy slightly differ. With four models applied in this study shown in Figs.
![]() | Fig. 5. The gNPT method determined the phase transition temperature of the TIP4P-2005 water model. Cyan region refers to 252 ± 5 K from the free energy calculations, yellow region refers to 249±3 K of the direct coexistence route.[21] |
![]() | Fig. 6. The gNPT method determines the phase transition temperature for the TIP4P-ICE water model. The melting point (272.2 K) proposed in previous work[30] is plotted as the dash line. |
We also give the results in two all-atom water models, as shown in Figs.
For the TIP4P-2005 water model, the obtained coexistence temperature is 250.0 ± 0.6 K, matching with Tc(DC) = 249 ± 3 K, Tc(FE) = 252±5 K, and Tc(DC) = 249.5±0.1 K from previous researches.[5,21] The temperatures in the three simulations with H = –4.25,–4.16,–4.06 come to almost the same value within 20 ns. The further 20 ns data are applied to estimate the coexistent temperature. The TIP4P-ICE water model is known to have a better performance in generating ice/water phase diagram. Its result Tc = 270.2±0.3 K, as shown in Fig.
To further explore the performance of this approach, we carry out simulations with different sizes of systems. Three gNPT simulations with 44800, 5971, 2030 water molecules are built, corresponding approximately to the sizes of 151 Å × 156 Å × 58 Å, 54 Å × 58 Å × 58 Å, 52 Å × 56 Å × 22 Å, respectively. These results are colored in red, as shown in Fig.
![]() | Table 1. The coexistence temperature of the mW, TIP4P/2005, TIP4P/ice water and LJ models as obtained from different methods (free energy calculations, Hamiltonian Gibbs–Duhem integration, direct coexistence technique, and the current GCE (gNPT) approach). . |
It has been found that the GCE (or gNPT) approach is very applicable of sampling metastable and unstable states of thermodynamics, as we did in this work sampling phase coexistence in the ice/water and the LJ solid/liquid systems. The GCE method can be applied to obtain directly the coexistent temperature of two phases even via a single individual simulation. Good precision of the approach can be reached even in a smaller simulation system, and it is better than (at least in comparable with) the other methods in the literatures with much less calculation time. Therefore, the GCE method is a highly efficient, accurate, and easy-to-use approach in determining the phase coexistent temperature.
[1] | |
[2] | |
[3] | |
[4] | |
[5] | |
[6] | |
[7] | |
[8] | |
[9] | |
[10] | |
[11] | |
[12] | |
[13] | |
[14] | |
[15] | |
[16] | |
[17] | |
[18] | |
[19] | |
[20] | |
[21] | |
[22] | |
[23] | |
[24] | |
[25] | |
[26] | |
[27] | |
[28] | |
[29] | |
[30] | |
[31] | |
[32] |